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Abstract. We consider a large class of two-lane driven diffusive systems in contact 
with reservoirs at their boundaries and develop a stability analysis as a method to 
derive the phase diagrams of such systems. We illustrate the method by deriving 
phase diagrams for the asymmetric exclusion process coupled to various second lanes: 
a diffusive lane; an asymmetric exclusion process with advection in the same direction 
as the first lane, and an asymmetric exclusion process with advection in the opposite 
direction. The competing currents on the two lanes naturally lead to a very rich 
phenomenology and we find a variety of phase diagrams. It is shown that the stability 
analysis is equivalent to an 'extremal current principle' for the total current in the two 
lanes. We also point to classes of models where both the stability analysis and the 
extremal current principle fail. 
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1. Introduction 

The physics of many non-equilibrium processes ranging from ionic conductors PQ and 
spin transport [2J to biological transport [3], [H [5J El [H E] can be described by driven 
diffusive systems. These are systems held out of equilibrium, for instance, by the 
environment driving a current of particles or energy through them. In the long-time 
limit, these systems can attain stationary states that are distinct from equilibrium states, 
and that escape description by conventional statistical mechanics. These systems have 
been widely studied for their fascinating behaviour, in particular, their propensity to 
exhibit boundary-induced phase transitions [91 [101 E] ■ 

Much work on driven diffusive systems has been carried out in one dimension — one- 
dimensional systems being of particular interest as they exhibit phase transitions with 
no equilibrium counterparts. A paradigmatic example is the totally asymmetric simple 
exclusion process (TASEP) whose phase diagram for the open boundary case has been 
studied in detail. Although in special cases such as the open boundary TASEP exact 
solutions are available [121 H31 EE], a useful approach for more general driven systems is 
to construct a mean-field theory. Remarkably, for the open boundary TASEP a simple 
mean-field approximation recovers the exact phase diagram j9j [15], HE] ■ Building on the 
mean-field theory, an "extremal current principle" has been proposed. This principle 
allows one to construct phase diagrams of one-dimensional open driven systems without 
explicitly solving the, generally non-linear, mean-field equations [TT | \T7 \ US]. 

More recently there has been much interest in extensions of the one- dimensional 
case to multiple lanes. A need for additional lanes arises naturally when describing 
many classes of systems. As an example, consider a biological transport system which 
describes the motion of molecular motors along a protofilament. Since molecular motors 
have a finite processivity they eventually detach and diffuse in the environment. To 
model the environment in a simplistic manner, a second lane in which the particles are 
diffusive and experience no exclusion can be introduced [19]. (For a different approach 
to tackle the same problem see [201 1211 (22] •) Two-lane models have also been used 
to describe the extraction of membrane tubes by molecular motors [23], macroscopic 
clustering phenomena [23] , spin transport [26] and various systems of oppositely moving 
particles [25l [271 EH EH]. In particular systems consisting of two coupled TASEPs have 
been the subject of several studies [321 [3U [351 EH1 EH [38] - More general multilane 
systems have also been studied [10~1 HI] and the hydrodynamics of coupled two-species 
systems has been considered in [311 HHl 1^2] - 

Given the broad applicability and interest in multilane driven diffusive models it is 
interesting to ask if there are generic principles, akin to the extremal current principle, 
which can be used to construct phase diagrams of such systems. In this paper we 
address this question for open boundary two-lane driven diffusive systems. We consider 
the broad class of models where particles can hop between lanes but where the motion 
within each lane is not influenced by other lanes. (For examples of the converse case 
in which particles cannot hop between lanes but where the motion within each lane 
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depends on the occupancy of the other lane, see [26, 30, 39J.) The dynamics in each of 
the lanes can have different characteristics, for example, asymmetric exclusion dynamics 
in one lane and diffusive dynamics in the other. We show that the phase diagram of such 
models can be constructed using a new stability analysis. The stability analysis is shown 
to be equivalent to an extremal current principle for two-lane systems and both methods 
may be used according to convenience. We illustrate the usefulness of the method for 
several systems and show that additional lanes generically lead to the appearance of new 
phases in the system. Finally, we argue that when the motion on each lane is influenced 
by the other lane the extremal current principle does not generically hold. 

The paper is organized as follows. In section we define the class of models that 
we shall study, and their mean-field dynamics. In section (|3| we present the stability 
analysis method for deriving phase diagrams of two-lane models. We illustrate this 
approach for a simple case in section (|4]) and then consider richer models in section ([5]). 
In section (|6]), we show an example of a more general class of systems where both 
our stability analysis and the extremal current principle may not apply. Finally, in 
section ^ we conclude. 



2. Model definition and mean-field equations 

2.1. The class of models 

Throughout the paper we consider a class of driven diffusive systems composed of two 
one-dimensional lattices (or 'lanes'). Each lane has L sites labelled i G [1, L] along which 
particles can hop. In addition, particles can hop between different lanes, the dynamics 
of which are thus effectively coupled. Each lane is connected at its ends to reservoirs 
which have specified particle densities. 

We will generally make the following assumptions: 

(i) The particles hop locally, i.e. particles can only hop from a site i to a neighbouring 
site in that lane, or to the same site i on another lane. (See for example figure [I]) 

(ii) The hopping rates of the particles within a lane depend only on the occupancies of 
their site and neighbouring sites within their lane (not on the occupancies of sites 
in the neighbouring lane). 

(iii) The average transverse flux of particles from one lane to the other increases with 
the occupancy of the departing lane and decreases with that of the arriving one. 

(iv) The reservoir densities are equilibrated, i.e. the densities are always chosen such 
that they produce zero transverse flux. 

(v) The dynamics in the bulk are translationally invariant. 
We note that ^ and (iii) are physically reasonable assumptions, whereas (jTTJ) , 



(iv) and Q are simplifying assumptions. In Appendix B.2 we discuss the effect of 



relaxing assumption (iv) whereas a full study of the effects of relaxing ^ (i.e. allowing 
direct coupling between the flows in different lanes) is beyond the scope of the present 
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paper. We do, however, argue in section (|6| that under generic conditions the methods 
developed in this paper may fail when this occurs and show a particular example where 
this happens. 

For a concrete example of a model in the class described above, see figure [TJ which 
illustrates a simple two-lane system. The dynamics of the bottom lane is that of a totally 
asymmetric simple exclusion process (TASEP). Namely, particles hop only to the right, 
with rate p, provided the arrival site is empty. In the upper lane, there is no exclusion 
and particles hop to the left and to the right with rates D~ and D + respectively, i.e. the 
particles perform biased diffusive motion. A particle on site i of the bottom lane can 
also hop to site i of the upper lane with rate d. A particle on site i of the upper lane can 
hop to site i on the lower lane with rate a provided that the arrival site is empty. Note 
that there are no explicit interactions between particles on different lanes — the hopping 
rate of the particles along each lane does not depend on the occupancies of the other 
lane. 

We refer to the model presented in figure [T] as Model I. We shall return to it 
throughout and use it as a template with which to illustrate our approach, before 
considering applications to other models towards the end of the paper. 
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Figure 1. Example of two coupled lattice gases in contact with reservoirs - Model I. 
The occupancies of sites i of the lower and upper lattices are denoted by t% and at, 
respectively. The bottom lattice has at most one particle per site and the particles 
can only jump to the right, at rate p (totally asymmetric exclusion dynamics). On the 
upper lattice, particles hop freely to the right with rate D + and to the left with rate 
D~ (biased diffusive dynamics). Particles on the bottom lattice hop with rate d to the 
upper lattice whereas those from the upper lattice hop with rate a to empty sites of 
the lower lattice. Both lattices are coupled to reservoirs with densities m, oi and r r , a r 
at their left and right ends, respectively. 



2.2. Mean- field equation 

To analyze the phase diagram of driven diffusive systems, it is often sufficient to consider 
a simple mean-field approximation, wherein one replaces the n-point correlation function 
by the product of n one-point correlation functions. For the average occupancies, one 
then ends up with equations of the form 

n = -j; + J l T _ 1 + K l ; cr^-jf + j^-Ki, (1) 

where T; and <7j are the average occupancies of the lower and upper lanes at site i\ J[ ', 
3° are the mean-field approximations of the current between sites i and i + 1 along the 
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lower and upper lanes, respectively. On the boundaries, the occupancies of both lanes 
equal those of the reservoirs, r h 07 and r r , a T (see figure 1). represents the mean-field 
value of the transverse currents between the two lanes. 
For Model I (see figure 1) the currents are given by 

Jl= m (l-T i+l ) 

J° = D+a, - D~a i+l 

Ki = -dr { + aai(l - Tj). (2) 

Although our approach could be carried out at the level of spatially discrete 
equations pj), it is convenient to take the continuum limit as detailed below. 

2.3. Continuum limit 

We take the continuum limit and keep terms up to second order in gradients. The 
mean-field equations are then given by 

t= -d x [J T -D T d x r} + K(r,a) (3) 
a = - d x [J a - D a d x a] - K(r, a). 

Here we have explicitly separated the mean field currents J T and J a into their advective 
parts, J T and J a , and their diffusive parts, —D T d x r and —D a d x o. The advective currents 
J T and J a are non-zero only when the dynamics on the corresponding lane are driven 
or equivalently when the hopping rates are asymmetric. Equations ^ have to be 
complemented with boundary conditions on the values of a and r which are imposed 
by the reservoirs: 

t(0) = t, r[L)=r r a(0) = a l a(L) = a r . (4) 



Following assumption (iv) we take the boundary densities to be equilibrated, that is 

K(T i ,a i ) = K(r r ,a r ) = 0. (5) 
We discuss the general case where ^ does not hold in Appendix B.2 



Our assumption (iii) about the hopping rates between the two lanes made in the 
previous section amounts to 

d T K < < d a K . (6) 

That is, the transverse current K increases and decreases with the occupancies of the 
departing and arriving sites, respectively. 

In the case of Model I (figure [l|, it is straightforward to deduce the currents from 
Eq.|2] 

J T =pr(l-r) (7) 

J a = VO 

K = ao~(l — t) — dr . 

Here v = D + — D~ , D T = p/2 and D a = (D + + D~)/2 and we set the lattice constant 
to unity. Thus, when D + = D~ , v = and the model reduces to that of [19]. Also note 
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that J T = J t (t) and J a = J CT (cr), which is consistent with assumption ([ii]), and that this 
choice of K satisfies condition (|6]). 

3. Stability analysis 

Constant (i.e. flat) profiles t(x) = r and a(x) = cr are solutions of the bulk 
equations ^ as long as they satisfy 

K(T ,a ) = 0. (8) 

We refer to such constant solutions as equilibrated plateaux. Because of our assumption 
(|6| on the dependence of K on r and a, the relation K(t , <to) = implies that o"o is an 
increasing function of r . Since the hopping rate within each lane does not depend on 
the occupancy of the other lane - see assumption (ii) - one can show that such solutions 
are dynamically stable. The details are left to Appendix A| where it is shown that 



perturbations about the plateaux 5t(x, t) = t(x, t) — r and 5a(x, t) = a(x, t) — a decay 
exponentially in time. On the other hand, we provide in section ^ an example of a more 
general two-lane system where assumption (ii) does not hold and where equilibrated 
plateaux are not necessarily dynamically stable. 

Let us now consider the steady-state version of the mean-field equations ^ 

d x (D T d x r) = d x J T (r)-K(r,a) (9) 
d x (D a d x a) = d x J a (a) + K(t, a) (10) 

along with the boundary conditions Q. Generally, constant profiles will not satisfy the 
boundary conditions Q on both sides since the left and right reservoirs may impose 
different densities. It is thus natural to ask whether one can construct steady-state 
profiles by matching two sets of plateaux that separately satisfy the left and right 
boundary conditions. As we now show, answering this question requires the study 
of the spatial stability of perturbations around the plateaux. We refer to perturbations 
as diverging if they grow as x increases and converging if they decrease as x increases 
and we refer to the corresponding plateaux as stable or unstable. It is important to 
note that from now on, any reference to stability or instability has to be understood 
with reference to the spatial dependence of a profile on x, since plateaux are always 
dynamically stable. 

In the example presented in figure [2j for instance, one needs the profiles to diverge 
away from the left plateaux and converge into the right ones as x increases. Since a is 
an increasing function of r in equilibrated plateaux, the perturbations for r and a away 
from the boundary plateaux must be of the same sign [|J 

To check that this indeed can happen we carry out a stability analysis of Eq. (]9 



and (10). Consider a small perturbation (e(x), r/(x)) around equilibrated plateaux so 



that t(x) = t + e(x) and a(x) = a + r](x) where r and do satisfy K(T ,ao) = 0. 
| Note that this assumes that the profiles connecting the plateaux are monotonous. 



Phase diagram of two-lane driven diffusive systems 



7 



n 



0~ r 



X 



L 



Figure 2. Plateaux imposed by the 'equilibrated' densities of left and right reservoirs 
are simple steady-state solutions of the mean-field equations at both ends of the system. 
It remains to be known whether the connecting shock is an acceptable steady-state 
solution. 



Expanding the steady-state mean-field equations to first order in e, rj gives 

D T d 2 x e = (d T J T )d x e - (d a K) V - (d T K)e (11) 
D a d 2 x V = (daJa)d x rj + (d a K)rj + (d T K)e. 

Here the derivatives with respect to r and a are taken at tq and ao respectively and 
similarly D T and D a are evaluated at r and a respectively Looking for solutions of 
this system of two second order linear differential equations of the form 

(e(x),r](x)) = exp(Xx)(e ,r] ), (12) 



equations (11) become 



, / X 2 D T - Xd T J T + d T K d a K \ ( e \ . .-, 

5 -d T K \ 2 D a - \d a J a - d a K ) Uo ) ■ ' ° 1 



For non-trivial solutions to exist, one needs the determinant of the matrix in (13) to 
vanish, that is 



A X (A) = (14) 



where 



X(A) = D T D a \ 3 - \ 2 {D T d a J a + D a d T J T ) + \(d a J a d T J T (15) 
+ D a d T K - D T d a K) + d a K d T J T - d T K d a J a . 

The spatial stability of the plateaux is then determined by the four solutions of the 
polynomial equation (14). In Appendix B we prove that the four solutions, which we 
denote Ao, Ai, A2, A3, are real and have the following properties: 

• There is a trivial solution Ao = which corresponds to shifting the values of the 
two plateaux while keeping K = 0, i.e. (e ,r/o) oc (—d a K,d T K); 

• There is exactly one additional solution such that eo^o > 0. The corresponding root 
Ai is the only root that may change sign as the parameters (r ,a ) are varied and 
its sign is given by the sign of x(0)- 
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There are two roots A2 < and A3 > with eo?7o < 0. 



• The root Ai is the only root of (14) relevant for the derivation of the phase diagram: 
it controls the convergence of perturbations away from constant densities. 

• The roots A2 and A3 only play a role for the boundary layers connecting the plateaux 
to non- equilibrated reservoirs (i.e. with K 7^ at the boundary) and not for the 
construction of the phase diagram (see Appendix B.2). 

The perturbations (e, rf) associated to Ai > thus grow as x increases (a diverging 
perturbation), while those associated to Ai < decrease (a converging perturbation). 
From equation ( |i~5| ) it readily follows that the change of spatial stability occurs when 

X(0) = d a K(r , a )d T J T (T ) - d T K(r , a )d a J a (a ) = 0. (16) 

Thus the perturbations (e, rj) grow or decrease with x when x(0) is negative or positive, 
respectively (see Appendix B.l). 

The method of constructing the phase diagram is then to determine the sign of 
x(0) for each set of equilibrated plateaux (imposed by the boundaries) and deduce 
whether the plateaux can be connected by appropriate perturbations. For example, the 
scenario presented in figure 2 requires a diverging perturbation from the left plateau and 
a converging perturbation at the right plateau. This can only take place if x(0) > for 
the left plateau densities 77, 07 and x(0) < for the right ones r r , ay. As we illustrate 
below, this line of reasoning provides a straightforward method to construct the phase 
diagram of multilane systems. 

Before turning to derive the phase diagram for several examples, let us show that 
x(0) is simply related to the derivative of the total current. In equilibrated plateaux, 
cr and r are related through K(t , cr ) = 0. The advective part of the total current, 
defined as 

Jtat = Jt + Ja 1 (17) 

may thus be expressed as a function of To whose total derivative is given by 

3 — ^tot(TO) °"o( r o)) — d T0 J t ot + d ao Jt t— — = d T0 J TQ — d ao J ao . 

Cir n <7Tn 0* A 



Equation ( 16 ) then implies 

X(0) = d Co K (19) 
dT 

Since under our assumptions (see Eq. rt6|) d ao K > 0, the sign of x(0) is the same as 
that of Therefore, the perturbation connecting different plateaux changes stability 
at an extremum of the advective part of the total current J tot . It is converging when 
< and diverging when > 0. 

dro dro 

Note that the case = requires a special treatment. Maximal current plateaux 
(where d ^2 ot < 0) are such that negative perturbations are diverging whereas positive 
perturbations are converging, so that profiles connected to these plateaux are always 
decreasing. Conversely, minimal current phases correspond to increasing profiles, i.e. 
t(x) and o~(x) increasing with x. 
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4. Construction of the phase diagram — a simple case 

To demonstrate how the stability analysis described above can be used to construct 
phase diagrams of different microscopic models we first study Model I when D + = D~ 
and thus v — (symmetric diffusion). Later we also consider the more general case 
where v ^ 0. 

Using equations ([TJ) and (19), we see that the stability of the equilibrated plateau 
is set by the sign of 

^=p(l-2r ). (20) 
dr 

This example is particularly simple because the advective part of the total current 



does not depend on o~o. Equation (20) implies that profiles with t < 1/2 are diverging 
(growing with x) whereas profiles with r > 1/2 are converging (decreasing with x). The 
case tq = 1/2 is marginal: negative perturbations are diverging whereas and positive 
perturbations are converging. This is consistent with our claim, in the previous section, 
that maximal profiles correspond to t(x) and o~(x) decreasing with x. The different 
profiles are illustrated in figure [3] 




Figure 3. Left: The left and right reservoirs impose boundary densities re and t t . 
Various possible steady profiles are presented, together with their stability. Right: 
Phase diagram constrained by the stability of the equilibrated plateaux presented on 
the left panel. 

The construction of the phase diagram is now straightforward. To do this we 
consider all possible matchings, four in total, between the different stability regions for 
left and right boundary plateaux. 

• When Tg > 1/2 and r r < 1/2, the only profile consistent with the stability analysis 
is flat at t(x) = 1/2 connecting to the boundaries at the ends of the system. This 
corresponds to the maximal current phase (MC phase). 

• When T£ < 1/2 and r r < 1/2, the only profile allowed by the stability analysis is 
a flat plateau with t(x) = T£ and a short boundary layer connecting to the right 
boundary density (LD phase). 
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• When n > 1/2 and r r > 1/2, similarly to the previous case, we have t(x) = r r and 
a short boundary layer connecting to the left boundary density (RD phase). 

• When T£ < 1/2 and r r > 1/2, the stability analysis is consistent with a shock 
profile connecting r(x) = T£ to r(x) = r T . The velocity of the shock, v s = 
{Jt{i~z) — Jr(Tr))/(Tr — Tg) is given by mass conservation. Therefore, the shock 
is steady for r r = 1 — Tt and is localized on the right or the left boundary for 
T r > 1 — ti and r r < 1 — re respectively. 

The corresponding phase diagram is drawn in figure [3] and is consistent with the 
results presented in [19]. This phase diagram reduces to that of the TASEP because of 
the lack of bias on the diffusive lane. In effect the diffusive lane adjusts its density to 
ensure that there is zero transverse flux between the lanes, therefore the phase diagram 
is determined by that of the driven lane i.e. the TASEP. 



5. Two-lane models with non-trivial phase diagrams 

In this section we consider more general two-lane systems which exhibit non-trivial phase 
diagrams. In these cases the phase diagrams are much richer than the TASEP coupled 
to a symmetric diffusive lane considered in the previous section, and we analyse in detail 
the phase diagrams of several such examples. 



5.1. TASEP and biased diffusion 

In this subsection we consider a TASEP coupled to a biased diffusive lane which 
corresponds to v > in Model I. Such a model is relevant to a biophysical situation 
where molecular motors diffuse asymmetrically after detaching from protofilaments, as 
for instance in models of cooperative extraction of membrane tubes [23J. 

Let us first show that in the presence of a bias in the diffusive lane, the phase 
diagram of Model I can become more complex and exhibit several new phases. As 
explained in section |3j the phase diagram can be deduced from dJ to t/dro. From the 
continuous mean-field equations ^ and one obtains 

J tot = pr(l - t) + va. (21) 

For equilibrated plateaux, K = imposes 

a 1 — t 

so that 

Jtot(ro) = — zr^— +pr (l - r ). (23) 
a 1 - r 

The values of tq where the stability changes are given by dJtot(7o)/dro = which reads 

(l-2r )(r -l) 2 = - — . (24) 

pa 
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Simple algebra shows that this equation always has an unphysical root at tq > 1. Fur- 
thermore, one can distinguish two regimes, according to the value of vd/pa, which exhibit 
distinct phase diagrams as we now discuss. 



Regime v > pa/ 2 Id: in this case there are no physical roots of equation (24). One 
finds that the plateaux are always diverging. Since all profiles are diverging, there is 
only one phase. One only observes a plateau with a density given by the density at the 
left boundary and a small boundary layer that connects to the right boundary. In this 
regime the advection in the diffusive lane is sufficiently strong that the left boundary 
always dominates. 



Regime pa / 27 d > v > 0: in this case there are two roots of (24), denoted tm and 
T m , which lie between and 1. Therefore the stability of the plateau changes twice: 
it is diverging for r close to and 1 and is converging in the intermediate region (see 
figure [4]). The smallest root tm is a local maximum of Jtot(To), so that profiles that are 
connected to the corresponding plateau are decreasing. On the contrary, the plateau at 
r m can only be connected to increasing profiles (see figure [4]). 

We now illustrate how the phase diagram can be constructed using figure |4j 
Once again, we look at possible matchings of the plateaux imposed by the boundary 
conditions. Both n and r r can be in the three different regions separated by tm an d T m . 
As indicated in figure [4j each of these nine cases has to be analyzed separately and the 
result is five distinct phases for any d,a ^ 0. Two of the phases disappear at d = 0. 

(i) Regions R(i), L(i) and L'(i): if the left and right densities are in the same stability 
region (a total of 3 different possibilities) there is one plateau set by one boundary 
and a small boundary layer that connects to the other boundary. For diverging 
regions (r r ,T£ < r M and T r ,r e > r m ) the left density controls the plateau whereas 
the right boundary dominates in the converging region (tm < T r , Tg < r m ). 

(ii) Regions mC(ii) and MCiii): when tm < T£ < r m and r r < tm, the stability 
implies a plateau at tm- Since this value corresponds to a local maximum of the 
current, we conclude that this region belongs to a maximal current phase. Similarly, 
when tm < re < r m and r r > r m , the stability implies a plateau at r m . This region 
therefore belongs to a minimal current phase. 

(iii) Regions L(iii) and R(iii): when T(. < tm and tm < r r < r m , the stability analysis 
predicts that the corresponding plateaux can be connected by an upward shock, 
from low to high density plateau. Since the total density p = r + a is locally 
conserved, the shock velocity is given by 

V s = [Jtot(Tr) - Jtatin)]/ (pr ~ Pi)- (25) 

Here, p r > pe and the plateau associated with the smallest current extends whereas 
the other recedes. The limiting case Jtotijr) = Jtotiji) defines a boundary line 
Tf (r r ) separating a 'right-dominated' phase R(iii) from a 'left-dominated' one L(iii), 
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corresponding to a discontinuous phase transition. On this line, a shock in the 
profile can in principle be observed. 

(iv) Regions L(iv) and mC(iv): when T£ < tm and r r > r m , the stability analysis 
predicts that a plateau of density t<? can be connected by a shock to a plateau of 
density r m which is then connected to the right density through a small boundary 
layer. This phase can be either 'left-dominated' or in a minimal current phase 
depending on the sign of the shock velocity. Since the shock is upward, the plateau 
with the lowest current dominates. A shock can in principle be observed on the 
transition line where Jtot( T e) — Jtotijm)- 

(v) Regions L(v) and R(v): when > r m and tm < r r < r m the stability analysis 
allows a downward shock to form in the system. Similar to the arguments regarding 
shocks presented above since pi > p r the plateau associated with the larger current 
controls the system. The phase transition line is given by Jtot(ir) = Jtat{Ti). 

(vi) Regions L(vi) and MC(vi): when > r m and r r < Tm the stability analysis 
allows a downward shock to form in the system from r(x) = Tg to a lower density 
given by tm- The latter is then connected to the right boundary with a short 
boundary layer. The scenario is depicted in figure |4| Again the larger current 
controls the system and there is a phase transition when Jtott/nu") = Jtotfa). 

Note that the contiguous regions L(i), L(iii), L{iv) are all controlled by the density 
at the left end of the system and therefore constitute a single, low-density, phase, with 
plateaux density smaller than tm- Similarly, the three contiguous phases L(v), L(vi) 
and L'(i) form a single high-density phase, in which the plateau density is larger than 
r m . Then, each of the groups of regions [R(i), R(iii), R(v)], [MC(ii),MC(vi)], [mC(iv), 
mC(ii)} defines a different phase. 

To summarize this subsection, we have constructed the phase diagram using a 
stability analysis. The resulting phase diagram has five distinct phases. In figure [4] the 
results of the analysis are compared to continuous time Monte-Carlo simulations and 
shown to agree very well. 



Note that we show in Appendix C that the above stability analysis leads to the 
same phase diagram as a generalization of the extremal current principle to the two- 
lane systems. For increasing profiles (ti < r r ) the steady plateaux are such that J t ot( r ) 
realizes its minimum for r G [t$,7v]. Conversely, for decreasing profiles (t# > r r ), the 
steady profiles are such that Jtot(T) realizes its maximum for r G [T r ,r^]. Which of the 
two methods one chooses to use is thus a matter of convenience. 

Next, we derive the phase diagram for two other two-lane models and compare the 
results to continuous time Monte-Carlo simulations. The derivation follows very similar 
lines to that presented above and we thus only outline the main differences. 

5.2. Two coupled TASEPs with advection in the same direction 

In this subsection we consider two coupled TASEPS with advection in the same direction 
but of different strength in the two lanes. 




Figure 4. Top-left: Current-density relation for the TASEP coupled to a biased 
diffusive lane. Top-right This figure represents a side view of steady density profiles 
that can be used to connect two reservoirs with densities r r and 77. According to the 
stability analysis, in the steady states, one can get diverging profiles for r < r M or 
t > r m and converging profiles otherwise. Note that the extremal densities have been 
spread out for clarity and their values are just illustrative. Bottom: Phase diagram 
of a TASEP coupled to a biased diffusive lane, for parameters p = 1, d = 0.01, a = 
1, D = 1, v = 1. The blue and red lines correspond to continuous and discontinuous 
phase transitions predicted by the theory while the black symbols were obtained from 
continuous time Monte-Carlo simulations. The grey dotted lines correspond to the 
local minimum and maximum of Jtot ■ Insets illustrate the stability discussion of the 
main text and the consequent profiles observed in the steady state. L and R stand 
for left- and right-dominated phase while mC and MC stand for minimal and maximal 
current phase. 
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For two coupled TASEPs with equilibrated boundaries, Harris and Stinchcombe 
have shown that if the extremal current principle applies, then the phase diagram 
should be the same as that of a single-lane TASEP with next-nearest neighbours 
interactions [35] . They backed up this assumption by simulating the profiles for several 
points in the phase diagram. Here we derive the phase diagram using the stability 
argument. Our results agree with those found in [35J and with continuous time Monte- 
Carlo simulations (see figure [6]). 

The model we study in this subsection is defined as follows. We consider two lanes, 
each consisting of a one-dimensional lattice with L sites. On the lower lane particles hop 
to a nearest neighbour empty site to their right with rate p. On the upper lane particles 
hop to a nearest neighbour empty site on the right with rate q. In addition particles can 
hop from site % on the upper lane to an empty site % on the lower lane with rate a and 
vice versa with rate d. To analyze the phase diagram we note that the currents along 
the two lanes are given by: 

J T =pr{l-r) ; J a = qa (1 - a) . (26) 

Because of the exclusion on both lanes, the transverse current K reads 

K = aa(l -r) - dr(l - a). (27) 

Using the K — condition, one can write the total current in an equilibrated plateaux 
as a function of To as 

da 



Ji 



tot 



To(l -r ) [p + q 



+ dr V 



28 



[a(l - r ) 

The total current as a function of r is shown in figure [5] for the most interesting case 
in which the curve has three extrema (depending on parameters there is also a less 
interesting case with only one extrema). This is the scenario we study below. 
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Figure 5. Left: Current density relation for two coupled TASEPs with a bias in the 
same direction. Right: Sketch of the stability of various steady-state profiles for two 
coupled TASEPs with a bias in the same direction. Note that the extremal densities 
have been spread out for clarity and their values are just illustrative. 
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As before the stability of the plateaux changes when dJ tot (r)/dr = 0. As can be 
inferred from figure [5] there are now three solutions which we denote, from smallest to 
largest, by tj^, r m and r| f . These lead to four distinct regimes of stability. When r < r\ t 
and r m < r < r| f the plateau are diverging and when t\ { < r < r m and r > the 
plateau are converging. Using a procedure identical to the one outlined above one can 
readily construct the phase diagram shown in figure [6j For example, consider phase L 
in the small T£ region. One then notes that for n < t\j and r r < t\j there is a plateau 
controlled by the left boundary with a small boundary layer on the right of the system. 
When Tg < and t]^ < r r < T m an upward shock connects the two boundaries. The 
system is then controlled by the smaller current which is given by J{r{) up to the line 
defined by J(r^) = J(r r ) and give rise again to a density controlled by the left boundary. 
Similarly when r e < r| f and r m < r r < r| f the scenario is identical but with the right 
plateau at r m and a small boundary layer at the right end of the system. Again it 
is easy to see that the system's density is controlled by Tg up to the line defined by 
J{ji) = J{i~m)- Finally, when n < rjj and r r > r| f the stability analysis allows for two 
upward shocks. One between ti and r m and one between r m and r r . It is straightforward 
to see that the plateau with r e dominates up to the line defined by J(r r ) = J(r^). 

The same arguments can easily be extended to construct the rest of the phase 
diagram shown in figure [6j The L (R) phase has a density controlled by the left (right) 
boundary and mC is a minimal current phase with density r m . The MC\ (MC2) phase 
is a maximal current phase with a density (r^). Finally, noting that Jir^) = Jir^) 
phase m corresponds to a shock phase with a density on the left of the shock and a 
density r\ f on its right. Note that there are two distinct phases controlled by the right 
boundary, one is a high-density phase in which r r is larger than whereas the other is 
intermediate between r]^ and r m . Similarly, there are two distinct phases controlled by 
the left boundary, one with a low density < and one with densities intermediate 
between r m and r^. 



5.3. Two coupled TASEPs with advection in opposite directions 

It is natural to look at the phase diagram of two TASEP lattices when the current in 
each of the lattices is directed in an opposite direction. Specifically we consider the 
two TASEP model considered above but where now particles on the upper lane hop to 
nearest neighbour empty sites to their left with rate q. 

Similarly to the case described above the currents on both lanes are given by 

J t = P t(1-t); J <T = -qa(l-a), (29) 



whereas the transverse current is still given by (27). The K — condition now yields 
for the total current J tot : 

Jtot = T {l-T ) (p-q— . ) ■ (30) 

V [a(l - T ) + dr \ 2 J 

In what follows we assume a > d (the regime a < d can be obtained by symmetry) and 
furthermore, to obtain a phase diagram which is significantly different from that of the 
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Figure 6. Phase diagram of two forward TASEPs with p = q = 1, d = 0.04, a = 1. 
The blue and red lines correspond to the continuous and discontinous phase transitions 
predicted by the stability analysis. The dotted grey line indicate the extrema of 
the current and the black dots are the results of the continuous time Monte-Carlo 
simulations. 




Figure 7. Current density relation, typical stability profiles for 2 TASEPs flowing in 
opposite directions for p = q = 1, d = 0.04 and a = 1. Note that the extremal densities 
have been spread out for clarity and their values are just illustrative. 



single lane TASEP, we take pa/d > q > pd/a. It is easy to check that when q is outside 
this range the current has only one extremum. The current in the regime of interest is 
shown in figure [7j As is evident there are now two points, denoted by tm and r m , where 
the stability changes. Plateaux are unstable for r < tm and r > r m and are stable for 
t m <r <r m . 

To construct the phase diagram we use the stability regimes (sketched in figure [7]) 
as discussed above. The results are shown in figure [8j Phases R and L correspond to 
phases which are controlled by the density of the right and left reservoirs respectively. 
MC corresponds to a maximal current phase with a density tm and the phase mC 
corresponds to a minimal current phase with a density r m . It is interesting to note 
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that inside phase R, say, as r r is increased, the total current in the system changes 
direction. This, however, happens in a smooth manner and is therefore not associated 
with a phase transition in the system. Note that as in the previous cases, there are two 
phases controlled by the left reservoirs, corresponding to low and high density phases 
(t£ < t m and T£ > r m , respectively). 

1.0 
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Figure 8. Phase diagram of 2 TASEPs with opposite advection directions, p = q = 1, 
d = 0.04 and a = 1. The area with a negative current total current is greyed. The 
dotted grey lines correspond to the extrema of the current while the black dots were 
obtained from continuous time Monte-Carlo simulations. 



6. A counter example — run-and-tumble model 

The derivations presented throughout the paper rely on assumption (ii) that the hopping 
rates of the particles within a lane depend only on the occupancies of their site and 
neighbouring sites within their lane. In this section we show that when the hopping rates 
on each lane depend on the occupancies of the others, implying interactions between 
lanes and a violation of assumption (ii), the situation is more subtle. In particular the 



results of Appendix A pertaining to the dynamical stability, do not always hold. To see 
this, below we explicitly look at a specific two-lane model with interactions. We show 
that in this model the equilibrated plateaux are not dynamically stable and therefore 
neither our stability analysis nor the extremal current principle hold. 

We consider a 'run and tumble' model with partial exclusion (see figure [9]). Particles 
at site i of the lower lane hop to the right with rate df while those on the upper lane 
hop to the left with rate d~ . To account for partial exclusion the rates df at which 
particles hop are given by 

df = v ± (l- T ^ + a ^ ) (31) 



P 

These rates decrease with the occupancy of the arrival site and are such that for any 
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site i the total occupancy Tj + o~i is always lower than p (which is treated as an external 
parameter). In addition, particles switch lane (or 'tumble') at rate a/2. 

The symmetric case without reservoirs was studied in [44] where it was shown than 
flat profiles are unstable at high enough density. The system undergoes phase-separation 
which results in an alternating profile with low and high density plateaux, separated by 
domain walls that perform random walks. The presence of advection (v + ^ v~) changes 
this picture slightly. At high enough density, equilibrated plateaux are again unstable 
and tends to phase separate. But as can be seen in figure [9j the domain walls propagate 
with non-zero velocity, according to a mass-conservation principle (25). 

Equilibrating the plateaux imposes r = cr and the current-density relation is thus 
very similar to that of the TASEP 

J(r ) = (*;+- t;>o (l - (32) 

The phase shown in figure [9] has however no counterpart in the TASEP phase diagram, 
which shows that the extremal current principle does not apply here. In addition, the 
current carried by the system can be shown to be unequal to the one imposed by the 
boundaries. 




50 100 150 200 

lattice site 



Figure 9. Left: Model of run-and-tumble particles on lattice. The interactions are 
set by the rates (31 1. Right: Result of continuous time Monte Carlo simulations for 
a = 1, v + = 12, v~ = 8, p = 100 and a system of L = 200 sites. Both left and right 
reservoirs tend to impose plateaux of density tq = uq = 26.5. At t = the system is 
prepared with plateaux tq = <tq = 26 but these plateaux are unstable and the system 
phase separates. Two successive snapshots of the total density Tj + Oi are plotted at 
t = 51 (black) and t = 71 (red). This shows the advection of the domain walls. 



7. Conclusion 

We have introduced a stability analysis which can be used to derive phase diagrams of 
driven diffusive systems. In particular, we have shown how the method can be applied 
to a class of two lane models and demonstrated its applicability using several examples. 
In some of them this allowed us to recover previously derived results with ease; we have 
also mapped out new phase diagrams for more elaborate models. Furthermore, the 
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method was shown to be equivalent to a extremal current principle and either method 
can be used according to convenience. 

The phase diagrams we have derived reveal some interesting physics. In the case of 
a totally asymmetric exclusion process (TASEP) coupled to a diffusive lane the phase 
diagram is unchanged from the single lane TASEP. However when a bias is introduced 
into the diffusive lane new phases emerge such as a minimal current phase and both low 
and high density phases controlled by the left boundary (figure 4). Similarly, when a 
TASEP is coupled to another TASEP a rich phase diagram including two left- and two 
right-boundary controlled phases, two maximal current phases and a minimal current 
phase emerges (figure 6). In addition for two coupled TASEPs with currents in opposite 
directions, a change in the direction of the total current may occur, although there is 
no associated phase transition (figure 8). 

It is important to note that we have also pointed out that the stability analysis, 
as well as the extremal current method, can fail for a broader class of models than 
outlined in this paper. For example, it can fail when the motion of particles on one lane 
is influenced by the motion on the other lane. To apply the method in such cases the 
dynamical stability of plateaux has to be considered in detail. 

Finally, we note that the method can be generalized to classes of multi-lane models. 
This will be discussed in detail in a future publication. 

Appendix A. Dynamical Stability of Equilibrated Plateaux 

In this appendix we show that for two-lane systems equilibrated plateaux are always 
dynamically stable provided the hopping rate along one lane does not depend on the 
occupancy of the other lane (a case 'without interactions' between the lanes). As shown 
in section [6] however, plateaux can be unstable if the aforementioned condition is not 
satisfied. 

Let us consider a perturbation around an equilibrated set of plateaux (r ,cr ): 

r(x) =r + J2 5r q e iqx - a(x) = a + ^ 5a q e iqx . (A.l) 

8 9 

At linear order in the perturbations, the Fourier modes decouple and the mean-field 
equations ^ yield 

d f5TA ~ M- M— ( " ^ + n . _ ) , (A.2) 



dt\SaJ \5aJ ' ' V -9rK -D a q 2 - iqj' a - d a K 

where J' T = d T J T and J' a = d a J a . The plateaux are dynamically stable if the eigenvalues 
of the matrix M have negative real parts. To see this, for simplicity, we first introduce 
a, b, c, d such that 

M= (" *) ■ (A.3) 



c d 

The eigenvalues A are then given by 



2A ± = {a + d)± v/(a - d) 2 + Abe . (A.4) 
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We denote by R(z) and I(z) the real and imaginary part of a complex number z, 
respectively Note that 

R(a + d) = -(D T + D a )q 2 - d a K + d T K < , (A.5) 

since all the terms are negative. For the real part of X^ 1 to be negative, it is thus sufficient 
to have 

\R(a + d)\> \R(V^) | with z = (a - d) 2 + 46c . (A.6) 
Using the relation 



R{Vz) = jmpA (at) 



we see that (A.6) is equivalent to 

2[R(a + d)] 2 > R[(a - d) 2 } + 46c + \z\ , (A.8) 
since be is real. Using this, a lengthy but straightforward algebra yields 

\a - d\ 2 - Abe + 8R(a)R(d) > \{a - d) 2 + 46c| . (A.9) 
Furthermore, since 

R(a)R(d) = D T D a q 4 -d T Kd (T K+q 2 (d a KD T -d T KD a ) > -d T Kd a K, (A.10) 
we have 8R(a)R(d) > 8bc and thus 

\a - d\ 2 - Abe + 8R(a)R(d) > \a - d\ 2 + Abc> \{a - d) 2 + 46c| , (A.ll) 

where the last inequality is the triangular inequality. Hence, the real parts of the two 
eigenvalues are negative and the plateau is dynamically stable. 

Appendix B. Eigenvalues and Eigenvectors of the Stationary Stability 
Problem 

In this appendix we present in more detail the stability of the perturbations around 
equilibrated plateaux and the role they play for the profiles. We first discuss in 



subsection Appendix B.l the solution of equation (14) 

=A X (A) (B.l) 

X (A) = D T D a X 3 - X 2 (D T d a J a + D a d T J T ) (B.2) 

+ X(d a J a d T J T + D a d T K - D T d a K) 

+ d a Kd T (J T )-d T Kd a (J a ), 

that controls the increase or decrease of the perturbations (e, rf) = e Xx (eo, 770) that satisfy 

n= f^D T -Xd T J T + d T K d a K \ fe Q \ /r> o\ 

{ -d T K \ 2 D a -\d a J a -d„K Uo J ' 1 1 



We then discuss in subsection Appendix B.2 the stability analysis near non-equilibrated 
reservoirs. 
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Appendix B.l. Eigenvalues 



First, one of the roots of (B.l ) is a trivial A = solution. The corresponding perturbation 
satisfies 



(e ,7?o) oc (-d a K,d T K). 



(B.4) 



Such perturbation corresponds to shifting the values of the two plateaux while keeping 
K = 0. Indeed, one can readily check that 



dK(r, a) = d T Ke + d a Kr] oc -d T Kd a K + d a Kd T K = 0. 



(B.5) 



We now show that there is only one solution (e , f]o) of equation ( |B.3[ ) that satisfies 

dws of 

^ ( DaVo 

D T e + D a i] D T i 

where 



£oVo > 0. To see this, we sum the two rows of equation (B.3) and divide by A to get 



D T e + D a r} 
d T J T 



D T e + D a r] ' 
d„J„ 



-B, 



(B.6) 



A =if and B =^7- (R7) 

If eo and 7]q have the same sign, A has to lie between A and B. One however notes that 
X (A) = (A - B)D a d T K; X (B) = (A - B)D T d a K. (B.8) 
Since d T K < and d a K > 0, one then has that 

• if A > B, x{A) < and X (B) > 

• if A < B, x{A) > and x{B) < 0, 

and thus x(min(74, B)) > and x( max (^5-B)) < 0. Since the coefficient of A 3 in x(A) 



is positive, there is only one root lying between A and B (see figure Bl). Since all 
this requires rjo and eo to be of the same sign, this means that there is only one root 
corresponding to such a perturbation, thus completing our proof. 

X(A) 

x[mm(A,B)] 

max(A, B) 




Figure Bl. Schematic plot of x(^)- There is at most one root between A and B. 



Note that for a marginal perturbation (A = 0) to exist one needs x(0) = 0, i.e. 

d a K(r , a )d T J T (T ) - d T K(r , a )d a J a (a ) = 0. (B.9) 
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This requires that A and B have opposite signs and thus the root that has vanished lies 
between A and B. This means that the only root that can change sign is associated 
with eo r] > C[§j 

We now consider the two other roots corresponding to eo^o < 0. They obey (see 
figure Bl) A2 < min(A, B) and A3 > max(A, B). If AB < then clearly A2 < and 
A3 > 0. In the case AB > we consider 

X(0) = D T Ad a K - D a Bd T K . (B.10) 



If A > and B > 0, then Ai > and A 3 > and, from (jBTOj), x(0) > 0. Since x(0) 
equals minus the product of the three roots of the cubic equation %(A) = 0, we deduce 
that A 2 < 0. Similarly, if A < and B < 0, we have A x < and A 2 < 0. Then x(0) < 
implies A3 > 0. Thus, in all cases A 2 < and A3 > and the sign of x(0) equals that of 
Al 



Appendix B.2. Connecting equilibrated plateaux to non- equilibrated reservoirs 



In this section we show that the perturbations that satisfy rjo^o < play a crucial role in 
the boundary layers that connect equilibrated plateaux to non-equilibrated reservoirs. 



Let us consider the situation presented in figure B2 The solution of the non-linear 
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Figure B2. Schematic representation of the system 



mean-field equations would lead to relations between plateau and reservoir densities 

r p = H 1 (r R ,a R ) ; a p = H 2 (r R , a R ) , (B.ll) 

where the R(P) superscript denotes densities at the reservoir (plateau). Since the 
plateaux densities are equilibrated, Hi and H 2 satisfy K(Hi,H 2 ) = K(p P ,p 2 ) = 0. 
Considering arbitrary perturbations 5t r and 5a R , the equilibrium condition (K = 0) 
imposes 

= 5p^[d^Kd^H 1 +d^Kd^H a ]+5 P ^[d^Kd^H 1 +d^Kd^H 2 ](B.12) 

§ Note that we have shown that the root associated to 770 eo > is the only one that can vanish. 
However, one has to be careful because the stability matrix is only diagonalisable at A 7^ so that the 
study of the extremal current phases has to be done separately. In this case, the degenerate A = 
eigenvalue is associated to a Jordan block. This is beyond the scope of this paper. 
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Since this has to hold for any perturbation, one gets 

d^Kd^H, + d pP Kd p nH 2 = (B.13) 

d p rKd pR H x + d pP Kd pR H 2 = (B.14) 

Last, this linear system of two equations for d p pK and d p pK has a non-trivial solution 
only if its determinant vanishes, which requires: 

d p RH 1 d p RH 2 - dpR.H2dpR.H1 = . (B.15) 

Next we construct a perturbation 8r , Sa R which is constrained to leave the plateau 
densities unchanged, namely Hi j2 (t r + 5t r , a R + 5a R ) = Hi j2 (t r , o~ R ). This requires 

= d^HxSr 11 + d^H^a 11 (B.16) 
= 3 t rH 2 5t r + d aR H 2 5a R (B.17) 

Again, this linear system admits a non-trivial solution iif the determinant vanishes so 
that 

d TR H 1 d tT RH 2 - d rR H 2 d aR H x = , (B.18) 



which holds due to Eq. (B.15). Finally, equations (B.16) and (B.17) show that the 
perturbation one has to look for is of the form 

{5t r , 5a R ) oc {-d aR H x , d^Hx) (B.19) 

or equivalently 

(5t r , 8a R ) oc (d a nH 2 , -d rR H 2 ) (B.20) 

This in particular means that 5t r 5o~ r < 0. Close to the left boundaries, such 
perturbations have to decay when x increases and they correspond to the eigenvector 
with components of opposite signs that has A < as discussed in the previous section. 
Conversely, the same reasoning for the right boundaries leads relevant perturbations 
with 5r R 5a R < and A > 0. This corresponds to the eigenvector with components of 
opposite sign and A > as discussed in the previous section. 

Concluding, for general systems with non-equilibrated reservoirs, there are 
boundary layers connecting the reservoirs to equilibrated bulk plateaux, whose behavior 
is dictated by the eigenvectors with e r7o < 0, one for each boundary. The connection 
between the two plateaux is then dictated by the sign of the remaining eigenvalue, that 
corresponds to e r] > 0. 



Appendix C. Connection with the Extremal Current Principle 

For single lane driven diffusive systems, studies based on domain wall theory have led 
to the formulation of the 'extremal current principle' [17]. In the multilane cases, it 
is also possible to show that for the transverse current K satisfying condition ^ and 
in the absence of interactions between the laneg|J[| our stability analysis shows that an 
'extremal current principle' applies. 

|| namely the hopping rate in one lane is independent of the occupancy on the other lane 
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Namely, for increasing profiles (ti < r r ) the steady profiles form an equilibrated 
set of plateaux such that Jtot( r ) realizes its minimum for r G [t^, T r \. Conversely, for 
decreasing profiles {ji > r r ), the steady profiles form an equilibrated set of plateaux such 
that Jtot(T) realizes its maximum for r G [r r ,r^]. The demonstration of both assertions 
being almost identical, we only consider here the increasing case. 

There are four cases to consider depending on whether T£ and r r correspond to 
stable or unstable plateaux, that is satisfy J t ' ot < or J' tot > (here quotes abbreviate 
total derivative with respect to r). Maximal current profiles are irrelevant since they 
correspond to decreasing profiles but minimal current profiles play an important role 
since they correspond to increasing profiles. When they exist, we denote by ... 
the densities for which J tot has local maxima. 



(i) J' tot {Ti) > and Jl ot {T r ) > 0. If J' tot {r) remains strictly positive on the interval 
[r^,r r ], Ti and r r are not separated by a minimal current phase and the steady-state 
profile corresponds to r(x) = Tg with a short boundary layer on the right. Since 
^tot( r ) > on [rf,r r ], the steady-state profile indeed minimize J tot on the interval 
[re,r r ). 

If J t ' ot vanishes on [n, r r ], the end-points are separated by at least one minimal 
current phase. A profile consistent with the stability analysis is given by a sequence 
of plateaux of densities in {77, r^} connected by up-going shocks. Mass conservation - 



equation (25 ) - shows that an up-going shock separating two density plateaux propagates 
towards the plateau with the largest current so that the plateaux with the smallest value 
of the current spread while the others recede. Since Jl ot (ri) > and J' tot {r r ) > 0, the 
smallest value among { Jtotiji), Jtotijm)} is indeed the global minimum of J to t in [rg, r r ]. 

(ii) J t ' ot (re) > and J' tot {r r ) < 0. is unstable and r r is stable. A profile consistent 
with the stability analysis is given by a sequence of plateaux of densities in {re, T % m , r r } 
connected by up-going shocks. Again, the plateaux with the smallest current spread 
while the others recede and the observed value of J to t is its global minimum in [re, r r ]. 

(Hi) J' tot {rg) < and Jt ot {r r ) < 0. and r r are stable. If there are no local minimum 
in [rg, r r ], the steady profile is a plateau at r = r r , J t ' ot remains strictly negative in [t£, r r ] 
and the current is indeed minimized over [t£, r r \. 

If J t ' ot changes sign, it has to do so at least twice and there are local minima 
r l m between rg and r r . A sequence of plateau of increasing densities r(x) G {T^,r r } 
connected by shocks is consistent with the stability analysis. Again, the plateaux with 
the lowest current spread while the others recede and the corresponding value of the 
current is the global minimum of J to t in [re, T r \. 

(iv) J' tot (Ti) < and J' tot {j r ) > 0. There is at least one minimum T l m between rg and r r . 
A sequence of plateau of increasing densities r(x) = T % m is consistent with the stability 
analysis. If there are several minima, the plateaux with the lowest value of the current 
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spread while the others recede. Since J' tot {T^) < and J' toXj {j r ) > 0, the global minimum 
of J tot is realized by one of the intermediate minimum and the value of J tot is thus indeed 
minimized on the whole interval. 

This concludes the list of 4 possible cases and shows that the current is always 
minimized between Ti and r r . Note that in all the above discussion, the minima of the 
current can be realized by several densities simultaneously. This leads to the coexistence 
of several plateaux with different densities. 

For decreasing profiles, the discussion is very similar, with two main differences. 
First minimal current profiles do not play any role while maximal current ones do, 
since only the latter are decreasing. Then, a downwards shock propagates towards the 
plateaux with the smallest current, so that the plateaux with the largest current spread. 
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